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Summary 


Hamilton's  principle  was  applied  to  derive  a  class  of  numerical  algorithms  for 
systems  of  ordinary  differential  equations  when  the  equations  are  derivable  from  a 
Lagrangian.  This  is  an  important  extension  into  the  time  domain  of  an  earlier  use  of 
Hamilton’s  principle  to  derive  algorithms  for  the  spatial  operators  in  Maxwell's 
equations.  In  that  work,  given  a  set  of  expansion  functions  for  spatial  dependences, 
the  Vlasov-Maxwell  equations  were  replaced  by  a  system  of  ordinary  differential 
equations  in  time;  but  the  question  of  solving  the  ordinary  differential  equations 
was  not  addressed.  Advantageous  properties  of  the  new  time-advance  algorithms 
were  identified  analytically  and  by  numerical  comparison  with  other  methods,  such 
as  Runge-Kutta  and  symplectic  algorithms.  This  approach  to  time  advance  can  be 
extended  to  include  partial  differential  equations  and  the  Vlasov-Maxwell 
equations.  Application  has  been  made  to  derive  a  second-order  accurate  algorithm 
for  the  linear  wave  equation;  the  dispersive  properties  of  the  algorithm  are  superior 
to  those  of  the  usual  second-order  accurate  explicit  or  implicit  algorithms. 
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Introduction 


In  the  description  of  plasmas  and  other  systems,  differential  and  integro- 
differential  equations  frequently  arise  which  are  derivable  from  a  Lagrangian. 
Examples  are  the  equations  of  motion  of  a  charged  particle  in  a  time-dependent 
electromagnetic  field;  the  Maxwell  equations;  the  Vlasov-Maxwell  system  of 
equations  for  a  collisionless  plasma;  and  wave  equations,  including  a  large  class  of 
nonlinear  wave  equations.  Other  examples  are  various  models  of  competition  or 
conflict,  such  as  the  time-dependent  Lotka-Volterra  equations  that  describe  a 
predator-prey  system  and  also  find  application  in  areas  such  as  mode  coupling  in 
plasma  physics  [1].  Thus,  development  of  algorithms  that  provide  reliable 
numerical  solution  of  Lagrangian  systems  of  equations  is  of  considerable  practical 
importance.  In  addition  to  requirements  of  accuracy,  numerical  stability,  and  speed, 
there  is  the  need  to  capture  the  important  qualitative  features  of  the  solutions,  for 
example,  the  phase-space  topology  of  solutions  of  particle  equations  of  motion.  In  a 
many-body  problem,  there  is  the  need  to  determine  average  quantities  accurately, 
even  though  most  details  may  be  unimportant.  For  example,  in  the  one¬ 
dimensional  Vlasov-Poisson  system,  a  quantity  of  particular  interest  is  the  energy  of 
the  electric  field,  which  is  an  average  quantity  determined  by  the  detailed  motion  of 
the  many  plasma  particles.  Usually,  numerical  algorithms  are  based  directly  on 
approximations  of  the  governing  equations.  However,  in  the  case  of  Lagrangian 
systems  of  equations  (or  their  Hamiltonian  equivalents),  Hamilton's  variational 
principle  can  be  used  as  the  basis  for  deriving  numerical  algorithms.  This  has  the 
potential  advantage  that  the  numerical  solutions  are  then  based  on  a  global 
principle  that  uniquely  determines  the  evolution  of  the  system.  It  is  in  contrast  to 
attempting  to  mimic  or  preserve  some  particular  properties  of  the  system  such  as, 
for  example,  conservation  of  the  Hamiltonian  for  an  autonomous  system,  or 
preservation  of  the  symplectic  nature  of  the  mapping  from  initial  data  to  final  data. 

Hamilton's  principle  was  introduced  for  deriving  approximation  schemes  for 
the  Vlasov-Maxwell  system  in  1968  [2].  Stated  briefly,  given  a  set  of  expansion 
functions  for  representing  the  spatial  dependence  of  the  electromagnetic  potentials, 
and  given  a  set  of  expansion  functions  for  representing  the  dependence  of  the 
particle  variables  on  the  phase-space  of  their  initial  conditions,  a  system  of  ordinary 
differential  equations  in  time  was  derived  whose  solution  determines  an 
approximation  to  the  solution  of  the  exact  Vlasov-Maxwell  system.  Specialization  to 
a  collection  of  point  particles  gave  particle-in-cell  simulation  schemes.  For  example, 
with  the  Vlasov-Poisson  system  in  two  cartesian  spatial  dimensions,  if  the 
expansion  functions  for  the  scalar  potential  are  linear  splines  ("tent  functions"), 
then  this  application  of  Hamilton’s  principle  leads  directly  to  the  standard  particle- 
in-cell  (PIC)  scheme  with  "area  weighting." 

It  is  well  known  that  time-advance  algorithms  are  very  relevant  to  accurate 
plasma  simulation.  Indeed,  Lewis,  Barnes  and  Melendez  [3]  raised  the  question  of 
whether  a  new  type  of  time-advance  algorithm  could  fundamentally  improve  the 
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quality  of  plasma  simulations.  The  original  application  of  Hamilton's  principle  to 
plasma  simulation  did  not  address  this  question;  no  time-advance  algorithms  were 
derived.  Nevertheless,  there  was  some  belief  that  practical  value  would  accrue  to 
having  time-advance  algorithms  based  on  Hamilton's  principle.  For  example, 
problems  with  numerical  stability  might  be  avoided,  or  it  might  be  possible  to 
reduce  the  effects  of  so-called  numerical  collisions.  Godfrey  derived  an  algorithm  to 
which  he  referred  in  an  unpublished  conference  proceedings  [4].  Although  he  found 
the  algorithm  advantageous  regarding  accuracy  and  stability,  he  also  found  it  too 
elaborate  and  costly  for  practical  use.  In  a  later  unpublished  conference  proceedings 
[5],  he  referred  to  a  practicable  approximate  implementation  of  the  algorithm  that 
was  sufficiently  accurate  for  a  two-dimensional  PIC  plasma  simulation  code  called 
CCUBE.  Details  of  the  algorithms  were  not  presented  in  these  proceedings  and  the 
work  appears  not  to  have  been  pursued  to  examine  the  merits  and  disadvantages  of 
time-advance  algorithms  based  on  Hamilton's  principle.  Eastwood  has  described  a 
"virtual-particle"  space-time  algorithm  for  electromagnetic  plasma  simulation  [6] 
whose  derivation  relies  on  Hamilton’s  principle.  However,  his  derivation  of  the 
time-advance  algorithm  is  inadequate  because  the  conditions  that  render 
Hamilton's  principle  equivalent  to  specifying  a  particular  solution  of  Hamilton's 
equations  are  not  addressed. 

Under  this  grant  (AFOSR  Grant  No.  F49620-92-J-0395),  an  initial  attempt  has 
been  made  to  determine  the  properties  of  time-advance  algorithms  based  on 
Hamilton's  principle  and  to  examine  their  merits  and  disadvantages  in  comparison 
to  other  approaches.  The  scope  of  the  activity  was  limited  to  ordinary  differential 
equations. 

The  Idea 


The  spirit  of  an  initial-value  problem  is  that  of  extrapolation:  Initial  data  are 
to  be  projected  to  a  later  time  according  to  the  action  of  governing  equations. 
However,  the  spirit  of  Hamilton's  priniciple  is  that  of  interpolation:  The  solution 
between  an  initial  time  and  a  final  time  is  to  be  determined  from  data  at  both  the 
initial  and  the  final  times  according  to  the  action  of  governing  equations.  The  key  to 
applying  Hamilton's  principle  to  derive  time-advance  algorithms  is  realizing  how 
to  bring  the  extrapolatory  nature  of  the  problem  into  consonance  with  the 
interpolatory  nature  of  the  principle.  There  are  two  steps  in  using  Hamilton's 
principle  to  derive  a  time-advance  algorithm.  First,  a  class  of  functions  with 
undetermined  parameters  is  chosen  for  approximating  the  solution  of  the  equations 
over  a  specific  time  interval.  For  example,  the  approximating  functions  could  be 
expansions  in  basis  functions  that  are  polynomials  in  time.  This  step,  in  one  form  or 
another,  is  part  of  any  derivation  of  a  time-advance  algorithm.  Second,  Hamilton's 
principle  is  applied  exactly  within  the  class  of  approximating  functions  to  obtain  the 
conditions  that  the  parameters  must  satisfy.  Those  conditions  involve  the 
unknown  value  of  the  solution  at  the  final  time.  However,  the  value  of  the 
solution  at  the  final  time  can  be  expressed  in  terms  of  the  initial  data  and  the 
undetermined  parameters.  Therefore,  the  conditions  on  the  undetermined 
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parameters  can  be  expressed  completely  in  terms  of  the  initial  data.  For  a  given  class 
of  functions  used  to  represent  approximate  solutions,  Hamilton's  principle  specifies 
a  unique  algorithm. 

Using  Hamilton's  principle  in  a  restricted  function  space  (the  space  of 
approximating  functions)  is  different  than  approximating  the  Lagrange  or  Hamilton 
equations  of  motion  directly.  If  the  function  space  were  not  restricted,  then  using 
Hamilton's  principle  would  be  equivalent  to  solving  the  exact  equations  of  motion. 
However,  if  the  function  space  be  restricted  to  a  set  of  approximating  functions  (for 
example,  polynomials  in  time  of  some  degree),  then  one  may  expect  the 
approximation  achieved  with  Hamilton's  principle  to  be  different  than  that 
obtained  by  applying  some  particular  approximate  method  to  the  equations  of 
motion  directly.  In  fact,  the  algorithms  based  on  Hamilton's  principle  are  equivalent 
to  a  Galerkin  approximation  of  the  differential  equations.  Hamilton's  principle 
determines  a  particular  weight  function  for  the  Galerkin  method  and  constrains  the 
function  space  in  which  Galerkin's  method  is  applied.  The  issue  of  which 
algorithms  are  best  is  complex,  and  it  can  only  be  resolved  by  detailed  numerical  and 
analytical  comparison  and  contrast  of  various  methods  with  one  another. 

Achievements 

The  following  progress  toward  developing  and  understanding  algorithms 
based  on  Hamilton's  principle  for  systems  of  ordinary  differential  equations  has 
been  achieved: 


•  A  general  formulation  has  been  developed  for  initial-value  problems. 
Extensions  of  the  formulation  can  be  made  to  a  wider  class  of  problems  (e.g.,  initial¬ 
boundary- value  problems). 

•  Algorithms  for  general  Lagrangians  have  been  defined  in  terms  of 
polynomial  approximating  functions. 


•  The  behavior  of  the  algorithms  for  small  values  of  the  timestep  has  been 
determined.  It  has  been  shown  that  the  algorithms  have  a  unique  solution  in  the 
limit  of  zero  timestep.  As  a  result,  they  can  be  implemented  in  a  numerically  stable 
manner  and  they  can  be  made  explicit  to  any  order  in  the  timestep. 


•  If  Hamilton's  principle  can  be  the  starting  point  for  deriving  algorithms 
that  are  optimal  in  a  useful  sense,  then  the  algorithms  should  have  analogs  of  every 
general  property  of  exact  Hamiltonian  and  Lagrangian  systems.  The  relation 


dH  dH  dL  .  . 

=  — —  IS  valid  for  ail  Hamiltonian  and  Lagrangian  systems.  For  that 
ut  at  at 


relation,  an  analog  has  indeed  been  demonstrated  for  computational  algorithms 
based  on  Hamilton's  principle  when  the  approximating  functions  for  the 
coordinates  and  momenta  are  polynomials  in  time. 
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•  Analytical  and  numerical  comparisons  of  algorithms  based  on  Hamilton's 
principle  have  been  carried  out  with  some  other  algorithms.  Among  the  algorithms 
used  for  comparison  were  symplectic  algorithms.  Symplectic  algorithms  mirror 
Hamilton's  equations  in  that  the  mapping  of  generalized  coordinates  and  their 
conjugate  momenta  from  one  time  to  a  timestep  later  is  canonical,  as  is  the  case  for 
the  exact  solution  of  Hamilton's  equations.  That  is  an  appealing  property,  although 
the  question  of  whether  it  is  the  best  property  to  strive  for  is  open.  After  all,  the 
main  objective  is  to  approximate  a  particular  solution  of  the  equations;  achieving  an 
approximation  that  represents  a  canonical  mapping  as  precisely  as  possible  may  or 
may  not  be  valuable  in  the  context  of  a  particular  problem.  It  might  have  been 
supposed  that  algorithms  based  on  Hamilton's  principle  would  also  be  symplectic; 
that  is  generally  not  the  case.  The  comparisons  of  symplectic  algorithms  with 
algorithms  based  on  Hamilton's  principle  indicated  that  the  properties  of  the  two 
kinds  of  algorithms  are  similar,  but  that  the  algorithms  based  on  Hamilton's 
principle  often  approximate  exact  solutions  better. 

•  Computer  code  using  Mathematica  and  C  has  been  written  to  enable 
systematic  numerical  comparisons  of  algorithms  based  on  Hamilton's  principle 
with  some  other  algorithms. 

In  addition  to  the  work  on  systems  of  ordinary  differential  equations, 
Hamilton's  priniciple  was  used  to  derive  a  second-order  accurate  algorithm  for  the 
linear  wave  equation.  It  was  determined  that  the  dispersive  properties  of  the 
algorithm  are  superior  to  those  of  the  usual  second-order  accurate  explicit  or 
implicit  algorithms.  The  stability  restriction  of  the  timestep  is  approximately  the 
same  as  with  the  usual  algorithms. 
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